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In patchy particle systems where there is competition between the self-assembly of finite clusters and liquid-vapour 
phase separation, reentrant phase behaviour is observed, with the system passing from a monomeric vapour phase to 
a region of liquid-vapour phase coexistence and then to a vapour phase of clusters as the temperature is decreased at 
constant density. Here, we present a classical statistical mechanical approach to the determination of the complete 
phase diagram of such a system. We model the system as a van der Waals fluid, but one where the monomers can 
assemble into monodisperse clusters that have no attractive interactions with any of the other species. The resulting 
phase diagrams show a clear region of reentrance. However, for the most physically reasonable parameter values of the 
model, this behaviour is restricted to a certain range of density, with phase separation still persisting at high densities. 



PACS numbers: 61.20.Gy, 64.75 .Yz, 64.60.De 

I. INTRODUCTION 

Self-assembly is a process by which a material, be it molec- 
ular, colloidal or macroscopic, spontaneously forms a well- 
defined structured aggregate. There are many examples of 
such self-assembling systems, be they biological (e.g. virus 
capsids, protein fibres and motors) or synthetic [e.g. surfac- 
tant micelles and supramolecular complexes ), although the 
latter are not yet able to replicate the exquisite control of self- 
assembly exhibited by biological systems. Indeed, some of 
the most remarkable synthetic self-assembling systems bor- 
row from biology; for example, the use of DNA and RNA to 
form complex nanostructures and devices. From the materi- 
als perspective, the advantages of self-assembly are clear, po- 
tentially enabling new materials and devices to be fabricated, 
particularly those with features on the nanoscale. 

For colloids and nanoparticles, there is a particular in- 
terest in developing particles that have 'patchy' anisotropic 
interactions, ■ with a view to extending the range of struc- 
tures that can be formed by self-assembly. Experimental 
methods for synthesising such particles, although in their 
early stages, are advancing rapidly. 1 1-21 Similarly, there is 
an increasing amount of computational work on patchy parti- 
cles, both as simple models for understanding biological self- 
assembly, 22-24 and to explore what might be possible with 
synthetic patchy particles and to learn the design principles 
for their successful self-assembly. 23 ^ 37 

One particular focus in the simulations of patchy par- 
ticles has been the formation of finite aggregates, be 
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they monodisperse ' 2 '' (similar to virus capsids) or 
polydisperse 33 ' 37-39 (micelle-like structures). One intrigu- 
ing result for these systems is that for certain param- 
eter ranges, the competition between self-assembly and 
liquid-vapour phase separation can lead to 'reentrant' phase 
behaviour. 29 ' 32 ' 33 ' 37 A phase transition is said to be reentrant if 
it involves the transformation of a system from one state into 
a macroscopically similar (or identical) state via at least two 
phase transitions through the variation of a single thermody- 
namic parameter (such as the temperature). 40 

Such reentrant behaviour is illustrated in Fig. 1 for a 
system that is designed to form monodisperse 12-particle 
icosahedra. The figure shows the average size of clusters 
as a function of the temperature and the angular width of 
the patches at the end of dynamical simulations that initially 
started from a monomeric gas. Values close to 12 indicate 
the formation of icosahedra, whereas very large cluster sizes 
indicate either a kinetic aggregate or the formation of a liq- 
uid droplet within a background vapour (and thus correspond 
to a liquid-vapour coexistence region). At very broad patch 
widths, the patches are relatively non-specific and the be- 
haviour tends to that of a Lennard-Jones fluid with phase 
separation between a vapour of monomers and a liquid of 
monomers. Conversely, at very narrow patch widths, the po- 
tential is so angularly dependent that the only thermodynam- 
ically stable phases are a vapour of monomers at higher tem- 
peratures and a vapour of clusters at lower temperatures (note 
that in these dynamical simulations, large aggregates form 
at low temperatures due to the slow kinetics, but this state 
is not a thermodynamically stable phase). At intermediate 
patch widths, we expect competition between self-assembly 
and phase separation to occur. For example, simulation results 
suggest that decreasing the temperature at a patch width of ap- 
proximately 0.5 radians and at a constant density involves two 
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FIG. 1. The mean size of clusters formed as a function of the reduced 
temperature and the patch width for patchy particles designed to form 
12-particle icosahedral clusters. These Monte Carlo simulations in- 
volved 120 particles simulated for 1.2 x 10 8 Monte Carlo steps at a 
constant number density of pa 3 — 0.15 and used the virtual move 
Monte Carlo algorithm. Results were averaged over 5 independent 
simulations. Decreasing the temperature at a patch width of approx- 
imately 0.5 radians involves two phase transitions; from a monomer 
gas through a liquid-gas coexistence region to a cluster gas. 



phase transitions: first going from a gas of monomers to a 
liquid-gas mixture and then to a gas of icosahedral clusters. 
This reentrance is driven by the lower energy of the cluster 
gas phase (due to the strong internal bonding within the icosa- 
hedra) compared to the liquid phase, which has a greater en- 
tropy. 

Clear evidence of reentrant phase behaviour has also been 
observed by Sciortino and co-workers in simulations of Janus 
particles with a single attractive patch. Certain parameter 
ranges lead to a competition between liquid-vapour phase sep- 
aration and micelle formation. • In the phase diagrams that 
they computed, the standard behaviour of an ever-larger two- 
phase region in the T-p phase diagram as the temperature is 
decreased below the critical point is perturbed by the forma- 
tion of micelles, shrinking the coexistence region and shifting 
it to higher density. However, what happens to this two-phase 
region at low temperature (does it end at a critical point, or 
persist to zero temperature?) is unclear, because the system in- 
stead formed a lamellar phase. In this paper, our aim is to ad- 
dress such questions by calculating complete phase diagrams 
for the fluid using a simple theoretical model of a system that 
can both self-assemble into clusters and exhibit liquid-vapour 
phase separation. 

Reentrant phase behaviour is of course not limited to patchy 
particle systems. Indeed, the classic example of such be- 
haviour is the nicotine-water mixture, which exhibits a closed 
loop in its temperature-composition phase diagram with both 
an upper and a lower critical solution temperature. Reducing 
the temperature whilst keeping the mole fraction of nicotine 
constant first leads to phase separation, and then to remix- 
ing driven by nicotine-water association. Liquid-liquid im- 
miscibility leading to such phase diagrams is often attributed 



to hydrogen bonding. • Complete mixing will occur if the 
temperature is sufficiently high so as to overcome the en- 
thalpic incompatibility of mixing. However, at sufficiently 
low temperatures, miscibility results from hydrogen bond- 
ing between different components in the system. Similar 
considerations have been used to explain reentrant phase be- 
haviour seen in polymer aggregation in water, in non-ionic 
micellar solutions and in nematic reentrant phases in liquid 
crystals. 47 



Although the systems of interest in this paper are not two- 
component mixtures, unlike most systems where reentrant 
phase behaviour has previously been observed, the behaviour 
can nevertheless in some ways be analogous to such sys- 
tems. The similarities arise due to the mixing of two states 
of the same component (monomers and associated clusters), 
rather than two components proper. An example of reentrant 
phase behaviour as a function of temperature in a system with 
a single effective density is ionic association. 48 ' 41 ' When di- 
electric constants are large, phase separation is driven by the 
usual enthalpic incompatibility. With low dielectric constants, 
coulombic-style phase separation between a low density (neu- 
tral cluster) structured gas phase and a high density (free ion) 
liquid phase occurs due to the long-range ionic interactions 
between the ionic components, rather than between the ions 
and the solvent. The formation of ionic clusters is somewhat 
similar to the cluster formation in the self-assembling patchy 
particle system that we focus upon, and so the reentrant be- 
haviour (in simple theories of these ionic systems, 'banana' - 
shaped coexistence curves are found ) may have similar ori- 
gins. 



In this work, we investigate a theoretical model which can 
be used to rationalise our earlier simulation results of self- 
assembling patchy colloidal particles and one which yields a 
phase diagram reminiscent of the work of Sciortino and co- 
workers. 33 ' 37 In particular, we are interested in constructing a 
phase diagram for single-component systems that can poten- 
tially form a vapour of monomers, a liquid of monomers and 
a gas of monodisperse clusters. Here, these clusters are for 
the most part chosen to be 12-particle icosahedra, but the 
approach can easily be generalised to other clusters, as in sub- 
section HID. We have chosen to model the chemical equilib- 
rium between clusters and monomers with explicit molecular 
partition functions, and to treat the vapour-liquid equilibrium 
as a familiar van der Waals fluid; this is a more straightfor- 
ward approach compared to using Wertheim theory, ■ which 
has already been used in the study of patchy particles. 28 We 
first look at just the monomer-cluster equilibrium in section II, 
and then determine the necessary changes to the pure van der 
Waals fluid model in order to be able to form the three-state 
'self-assembling van der Waals fluid' and calculate appropri- 
ate phase diagrams for this system in section III. 
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II. MONOMER-CLUSTER EQUILIBRIUM 
A. Partition functions 

We first construct the molecular partition function of an 
icosahedral cluster in a manner similar to that done previously 
for a simple model of tetrameric protein complexes. 23 ' 52 

The total molecular partition function of the icosahedral 
cluster is given by 



<7dus = <7tr<Zrot9vib exp [-EjkT] 



(1) 



where E c is the ground state energy of the cluster. As we are 
comparing results to classical simulations, and since we are 
modelling colloidal particles, it is reasonable to use classical 
expressions for the various partition functions. The transla- 
tional partition function is given by 
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where A is the de Broglie thermal wavelength, the mass of the 
icosahedron is 12m (with m being the mass of a monomer), 
and N mon and iV c i us refer to the number of monomers and the 
number of clusters in the system, respectively. For the vol- 
ume term, we choose to use the free volume approximation, 
which is one generalisation of the van der Waals expression 
to multiple components. The volume available to clusters is 
thus V — BN c i us — bN mon , taking into account the reduction 
in the total volume due to the presence of other clusters and 
monomers. Here, B and b are free volume coefficients (rather 
than hard-sphere second virial coefficients) for clusters and 
monomers, respectively. The choice of values we assign to B 
and b is discussed in subsection II B. 

The moment of inertia about the principal axes of an icosa- 
hedron is given by I c = 80toct 2 /(5 — V5), where m is the 
monomer mass and a is the diameter of the monomers. The 
symmetry number of an icosahedral cluster is 60. We can 
therefore calculate the rotational partition function as 
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There will be (6x12 — 6) overall vibrational modes of the 
cluster. For simplicity, we assume that (3x12 — 6) probe the 
radial part of the potential, and 3 x 12 probe the orientational 
degrees of freedom, and that all the modes of each type have 
the same spring constant. Hence, we can write the vibrational 
partition function as 
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To obtain expressions for uj m a and w ang , we need to assume 
some form for the interparticle potential. We choose the po- 
tential that we used in our simulations of the self-assembly of 
icosahedra. • The potential is based on the Lennard-Jones 
form, but where the attractive region is modulated by an an- 
gular term that measures how well patches point at each other, 
and is given by 
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where ry is the inter-particle vector connecting the centres 
of the two particles i and j, n* is its magnitude, i~2j and flj 
are the orientations of the particles i and j, the Lennard-Jones 
potential is 
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where cr pw is the patch width (we measure er pw in radians), Okij 
is the angle between the patch vector of patch k on particle 
i and the interparticle vector r,,-, and fc m ; n is the patch that 
minimises the magnitude of this angle. Hence, two particles 
interact only through a single pair of patches. 

We assume that £j ra( j and cj ang are simply related to the sec- 
ond derivative of the radial and angular parts of the above 
anisotropic modified Lennard-Jones potential. Namely, they 
are given by 
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where 7 mon is the moment of inertia of the monomer (given by 
J mon = (2/5)mcr 2 for a spherical monomer), and 
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Although these approximations appear somewhat crude, this 
level of description is reasonable given some of the other ap- 
proximations used elsewhere. Furthermore, using the exact 
form of the vibrational partition function would be unlikely to 
have a significant effect on the overall phase behaviour. 

Finally, the Boltzmann factor in Eq. 1 reflects the addi- 
tional energy obtained as a result of clustering - that is to say, 
the energy of bonding between the patches of the monomers 
when arranged in an icosahedron. Ignoring next-neighbour 
interactions, the ground state energy of the icosahedron is 
E c = — 30e, since each of the 12 monomers in an icosahe- 
dron is bonded to five other monomers with bond energy e. 
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Overall, the icosahedral cluster molecular partition function 
is given by 
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The monomers have both rotational and translational de- 
grees of freedom, but not generic attractions at this stage. 
Therefore, the molecular partition function of the monomer 
is given by q mon = q a> mon x q roL mon , where 
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Note that the symmetry number s of a monomer is 5, since 
there are five identical patches on each monomer, and the 
same volume term occurs as in the cluster translational 
partition function because we are using the free volume 
approximation. 48 The overall monomer molecular partition 
function is therefore given by 
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B. Estimation of free volume coefficients 



We choose to use the free volume approximation sum- 
mation to take into account the effects of excluded volume. 
Although there are more accurate functions which describe 
such effects better, we deliberately wish to keep the par- 
tition functions as simple as possible, as we are interested 
in the generic properties and phase behaviour resulting from 
the competition between self-assembly and phase separation, 
rather than the detailed behaviour of a specific system. More- 
over, we want a form consistent with the van der Waals de- 
scription of the fluid, and one that works reasonably well both 
at high and, to a degree, at low system densities. 

Within the free volume approximation, the maximum 
density in the system will equal the reciprocal of the free vol- 
ume coefficient b, p max = 1/b. Throughout this paper, all den- 
sities p refer to number densities, p — N/V, rather than mass 
densities. We choose to set the maximum density achievable 
in the system to 1 er~ 3 , and hence b = a 3 . This choice is 
reasonable both in terms of the comparison of the maximum 
packing fraction with that for random close packing, and of 
the position of the van der Waals critical point with that of the 
Lennard- Jones fluid. 

The relative values of B and b should be proportional to 
the cluster and monomer sizes. To determine the size of the 
cluster relative to that of the monomer, we assume that clusters 



are spherical objects, but such that we take their icosahedral 
character into account. The circu mradius of a n icosahedron 
with edge length a is r = ((1/4)^10 + 2^5 ~ 0.95 a; if the 
icosahedron is built of monomelic spheres of diameter er, then 
a = a. However, we also need to take into account the fact 
that these monomers protrude out of the icosahedron itself, 
and so another hard-core radius of these constituent spheres, 
<r/2, must be added to the circumradius of the icosahedron to 
estimate the overall cluster radius R, giving R = 1.45 a. 

The volume occupied by this sphere relative to the volume 
occupied by a single particle is R 3 /(a/2) 3 = 24.4. This is 
therefore a suitable measure of how much bigger the icosahe- 
dral cluster is compared to a single patchy particle, and sug- 
gests that an appropriate ratio of cluster and monomer coef- 
ficients is B/b w 24. We consider this case most compre- 
hensively, and also a few other values in order to characterise 
more fully the types of phase diagram that may arise. This ra- 
tio of B jb is specific to the icosahedral system described, and 
would take values different from the ones estimated here if the 
system were allowed, for example, to adopt a different cluster 
shape (for the same number of monomers) or to assemble into 
a different oligomeric state. 



C. Clustering transition 

We assume that the only relevant reaction in the system 
is 12A =± A12, which is to say that there are no partially 
constructed clusters whatsoever. We justify this assumption 
from simulations, which show that the free energy of a 12- 
particle cluster is considerably lower than any intermediate 
size bigger than a monomer. At equilibrium, the chemical 
potentials of the individual species, appropriately weighted by 
their stoichiometric coefficients, will be equal; in this case, 
12^ mon = /i c i us , where pi = —kT\n(qi/Ni). We choose to 
define the centre of this equilibrium as the point at which the 
probability of a particle being in a cluster or a monomer is 
equal, which imposes the condition p mon = 12p c \^ s , where 
Anon = N mm /V and p dm = A clus /F. We denote the temper- 
ature at which both conditions are fulfilled as the clustering 
temperature, T C ] US . 

We are not able to solve this system of equations analyti- 
cally, and so we use numerical minimisation of the difference 
in chemical potential using the density constraint to obtain 
T c ius as a function of density, as plotted in Fig. 2. The trans- 
formation between the two states is continuous (see Fig. 3) 
and the T c \ us (p) line is not a phase coexistence curve, but sim- 
ply the centre of a chemical equilibrium that indicates where 
clusters and monomers are equally probable states. 

We first consider the ideal limit (6 = and B = Q) and 
see from Fig. 2 that at low temperatures, the cluster gas is 
more stable than the monomer gas. Although the clusters 
have a lower entropy as a consequence of the lower number 
of translational degrees of freedom compared to the monomer 
gas, at sufficiently low temperatures, the clusters will have a 
lower free energy, because their lower entropy is more than 
overcome by their lower energy due to bonding between the 
patches. Also, the equilibrium favours the cluster state as the 
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FIG. 2. Tdus (p) for the transition from monomers to icosahedral clus- 
ters for different ratios of the free volume coefficients, B/b. For the 
non-ideal systems, b/a s = 1, and dotted lines show plSl^a 3 , the 
maximum density at which particles can have an equal probability of 
being monomeric or in a cluster. <j pw = 0.5. 



density increases (i.e. T c i us increases with p), since the effec- 
tive number of particles in the gas phase is lower. 

Although the above are the most important basic features 
of the monomer-cluster equilibrium, there are a few more fea- 
tures that we should consider. For example, the model predicts 
that at high temperatures, the cluster gas becomes favoured. 
This behaviour is an unphysical consequence of applying the 
harmonic approximation in a temperature range where one ex- 
pects all bonds to be dissociated. However, this deficiency in 
the model is not a problem for the current study, as we are 
interested in temperatures that are well below where this phe- 
nomenon occurs. 

Figure 2 also shows the effect of introducing excluded vol- 
ume interactions on the monomer-cluster equilibrium. There 
is no maximum density in principle for the ideal system 
(b = 0, B = 0), and clusters are always stable to the right of 
the clustering transition line. The first effect of the excluded 
volume is to limit the maximum density. For the B/b = 24 
case, the maximum density of monomers is p mon , maxC 3 = 1, 
and the maximum particle density associated with the clusters 
is 12p c i US; maxC 3 = 1/2. For the B/b — 12 case, both the 
monomers and the clusters have a maximum particle density 
of 1. 

At low densities, the excluded volume interactions have lit- 
tle effect and the T c ius(p) lines initially follow that for the 
ideal case. However, as the density increases, the lines are 
displaced upwards with respect to the ideal case because the 
volume term V — BN C \ US — bN mon in the translational partition 
functions of both monomers and clusters decreases, and this 
destabilises the monomeric state more because it has a greater 
number of translational degrees of freedom. 

A typical dependence of the fluid composition on the den- 
sity is shown in Fig. 3 for a system in which B/b > 12. The 
fraction of particles in clusters increases rapidly as the system 
passes through the clustering transition on increasing the den- 
sity, and as p = 12/9 c i usmax is approached, this fraction is usu- 
ally close to unity. However, when clusters exclude more vol- 
ume than 12 monomers (which is the usual physical situation), 




FIG. 3. The relative proportion of monomers and clusters in the equi- 
librium mixture is shown as a function of the density at kT/e = 0.2 
for B/b = 24, <r pw = 0.5. The function (12p c i us — Pmm)/p takes 
the value —1 when only monomers are present, +1 when only clus- 
ters are in the system, and at the centre of the clustering transition 
(i.e. on the T c i us (p) line). Vertical dotted lines indicate the maximum 
particle density associated with clusters 12p c ius, max, the maximum 
density at which the monomers and clusters can have an equal proba- 
bility pm'ax and the maximum monomer density p raoll _ max , as labelled. 
The solid line is calculated from the partition functions, whilst the 
dashed line is the extrapolation in the regime where the available 
volume is held at zero. 



i.e. when B > 12b, achieving densities beyond 12p c i usmax 
is only possible if the monomeric state again becomes popu- 
lated, because monomers pack more efficiently (for example, 
when B/b = 24, 12 monomers take up only half the vol- 
ume of a single icosahedral cluster). Hence we see in Fig. 3 
that as the density increases further, the fraction of monomers 
increases, becoming the dominant species beyond pm ax , the 
maximum density for which particles can have an equal prob- 
ability of being monomeric or in a cluster. In Fig. 2, the 
dotted vertical lines correspond to p^lx, and so for B/b = 24, 
monomers dominate to the right of the line. We also note that 
in practice, it becomes numerically extremely difficult to cal- 
culate the equilibrium densities of monomers and clusters be- 
yond p = 12p c i USi max , because the V — BN C \ US — bN mon term in 
the translational partition function becomes very close to zero 
and the chemical potentials of the clusters and the monomers 
both become extremely large. However, it is relatively easy to 
extrapolate beyond 12p c i us>max using the approximation that 
the available volume V — BN c \ as — bN mon is zero, as illus- 
trated in Fig. 3. 

We note here that when considering the combined self- 
assembling van der Waals fluid in section III, for the most part 
we need not worry about some of these phenomena that occur 
when p > 12p c i us max , because phase separation intervenes 
before this density is reached. 
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III. SELF-ASSEMBLING VAN DER WAALS FLUID 

Having now developed an approach to calculate the equi- 
librium between a vapour of monomers and clusters, we now 
wish to calculate the phase diagram of a 'self-assembling' van 
der Waals fluid, where monomers can either assemble into 
inert clusters (with no attractive interactions) or condense to 
form a liquid. We have chosen to use a van der Waals fluid 
as a starting point because, although we could use a more so- 
phisticated description of the monomer and liquid states, we 
are mainly interested in the underlying fundamental physical 
behaviour of the system rather than in trying to describe the 
behaviour of any one particular experimental system. 



A. Partition functions 

To model clusters, we use the cluster molecular partition 
function from Eq. (10) above. The canonical partition func- 
tion for a van der Waals (vdW) monatomic fluid is given by 



1 /2nmkT\ {3N/2) 



aN' 



2 1 



VkT 



, (14) 



where a and 6 are positive constants. This expression is anal- 
ogous to that for an ideal gas, but with an excluded volume 
term (Nb) and an attractive exponential term. To model the 
vapour and the liquid, we use the molecular partition func- 
tion for a van der Waals fluid derived from Eq. (14), except 
for two modifications. First, we include an additional term 
in the volume expression in order to take into account the 
exclusion of monomers by clusters as well as by monomers 
themselves. Using the free volume approximation leads 
to the same available volume term as in Eq. (10), namely 
V—BN a i m — bN mon . Furthermore, we multiply this molecular 
partition function by q Iol of the same form as in Eq. ( 1 2) to take 
into account the fact that the monomers are now particles with 
orientational degrees of freedom. There are no attractive inter- 
actions between clusters, or between clusters and monomers, 
to reflect the fact that in our patchy particle model, no patches 
are exposed on the surface of the clusters (we ignore partially 
formed clusters); all patches are instead used in the internal 
bonding. 



B. Coexistence curve 

The coexistence curve, i.e. the curve where the difference 
in the Helmholtz energy between two phases is zero, is often 
calculated using a Maxwell construction, whereby a straight 
line is constructed in the p-V phase diagram over the me- 
chanically unstable region such that the areas of the regions 
defined by the van der Waals loop and this line are equal. 
However, in many interesting cases, van der Waals loops do 
not exist, 48 ' 49 and it would be useful to have a more generic 
(and, indeed, easier) approach. We can calculate the chemical 
potential of species j using fXj = (dA/dNj)y,T.N', where 



N' refers to all components other than j, and the pressure as 
p = —(dA/dV)T,N- Plotting the chemical potential against 
the pressure parametrically allows us to determine binodal 
points where the curves cross. Whenever the chemical po- 
tential defines two (or more) different pressures, the stable 
phase invariably has the largest pressure, which maximises 
the entropy and is consequently favoured by the second law 
of thermodynamics. 48 ' 57 

As in section II, we assume that the only relevant reac- 
tion in the system is 12A ^ A%2. We then seek to find 
the coexistence curves for the system where van der Waals 
monomers and clusters are in equilibrium. This is achieved 
by setting fi c \ us = 12/i mon and constructing p-fi curves to cal- 
culate binodals. Since the density of clusters is a function 
of the density of van der Waals monomers, and vice versa, 
we use numerical techniques for solving the system of equa- 
tions. To find p as a function of /i, we impose a chemical 
potential fi and numerically minimise both |/i c i us — 12/^ mon | 
and |/icius/12 + jUmon — 2/i| as functions of p mon and p cms , all 
under the constraints that densities are positive and that the 
available volume (V — BN c i us — 6A c i us ) is non-negative. The 
first minimisation ensures that the system is at equilibrium 
(Mcius ~ 12/i m on), whereas the second fixes the individual 
chemical potentials to the imposed value (/i c i us ~ 12/i and 
A*mon ~ m)- We only admit answers which have a combined 
difference equalling zero to within a small error term of 10~ 6 , 
although in practice the error term is generally several orders 
of magnitude smaller. Having now calculated p c i us and p mon , 
we can calculate the pressure. Example p-v and p-fi curves 
are given in Fig. 4. By performing this analysis at a large 
number of relevant temperatures, the full phase diagram can 
be calculated. 



C. Stability criteria 

To include spinodals on the phase diagram, we can bracket 
regions of mechanical stability by finding where the gradi- 
ents of the p-v curves change sign. Notice from Fig. 4(c) that 
there are parameter ranges where there are two regions of me- 
chanical instability (AB and CD), separated by a region of 
mechanical stability (BC). In this system, however, we also 
have to check for compositional stability, i.e. whether fluctu- 
ations away from the equilibrium composition of monomers 
and clusters are stable. To derive a condition for composi- 
tional stability, we start from the fundamental equation for the 
Helmholtz energy, which may be written as 

dA= - p dV - S dT + /Union dA m0 n + Mdus dN dm . (15) 

However, the number of monomers and clusters is not inde- 
pendent, but N = 12A dus + A mon and dA = 12 dA dus + 
dA mon . Hence, we may rewrite the above equation as 

dA= -pdV - SdT + p mm dN + 



+ (Mclus — 12/i m on) dA dus . 



It follows that 



OA 



dN c] 



Mcius — 12/i 



(16) 



(17) 



clus / N V T 



7 




FIG. 4. The p-v ((a), (c) and (e)) and p-fi ((b), (d) and (f)) curves for the self-assembling van der Waals fluid at reduced temperatures of 
(a) and (b) kT/e = 0.21, (c) and (d) kT/e = 0.17 and (e) and (f) kT/e = 0.02. Regions of mechanical stability are coloured blue, those 
of mechanical instability red, and the mechanically stable, but compositionally unstable, region of back-bending is coloured green. Binodal 
points are marked by violet asterisks connected by a tie line. Spinodal points, where the derivative of the pressure with respect to the volume 
is either zero or infinity, are also shown with the labelling and colour-coding matching that used in Fig. 5. v = V/N, a/ea 3 = 1, b/a 3 = 1, 
<7 pw = 0.5,B/6 = 24. 



which is of course zero at equilibrium. The condition for com- 
positional stability is hence that 



d 2 A 



( <9(/Jrius - Virion) 

V cW clus 



> 0. 



(18) 



N, V,T 



This stability criterion can be evaluated at every equilibrium 
point of the combined system. The system is compositionally 



unstable exactly along the line BC. 



D. Phase diagrams 



We are most interested in computing phase diagrams for our 
system in the region of parameter space where reentrant be- 



8 




FIG. 5. The T-p phase diagram for the self-assembling van der Waals fluid for B/b = 24. The thick lines give the binodals for the system. 
Also represented are lines which correspond to points A, B, D and E in the p-v plots in Fig. 4, where the derivative of the pressure with respect 
to the volume is zero or ±oo, where they are true limits of stability. The binodals and spinodals for the standard, pure van der Waals fluid 
are plotted for comparison (dotted lines). Finally, the clustering temperature (T c i us ) is also plotted (thick dotted line), a/ea 3 = 1, b/a 3 = 1, 
<7 pw = 0.5. 



haviour might be expected to occur. We choose a value for the 
patch width from the region where simulations suggest reen- 
trant phase behaviour occurs (Fig. 1), namely er pw = 0.5. We 
choose the free volume coefficients B and b in line with sub- 
section IIB, such that b/a 3 = 1, and use the case B/b = 24 
as our main example. We then choose the van der Waals pa- 
rameter a such that the critical van der Waals temperature 
kT c = 8a /27b is somewhat greater than the clustering tem- 
perature at the critical density. This is true for our choice 
a/ea 3 = 1, for which fcT cl -; t /e ~ 0.30 for the pure van der 
Waals fluid and kT c i us (p cr i t ) / e sw 0.23 in the absence of liquid 
formation (i.e. a = 0). 

As in subsection II C, we can also calculate the clustering 
transition temperature, T i us (p), although we should note that 
where T c i us (p) lies in the two-phase region, phase separation 
of course occurs instead of clustering. With liquid forma- 
tion being favourable at high densities, however, the clustering 
transition curve does not simply increase in temperature as the 
density increases, but beyond pa 3 s» 0.15 instead gradually 
decreases, such that it reaches p^lx at kT/e — 0, as shown in 
Fig. 5. 

The full phase diagram for B/b = 24 is depicted in Fig. 5. 
Reentrant phase behaviour is readily observed in this T-p 
plane for 0.05 J$ pa 3 ;$ 0.5. As the temperature decreases at 
constant density within this range, the system first undergoes 
a transition from a monomeric fluid to a two-phase liquid- 
vapour coexistence region, and then as the temperature is fur- 
ther decreased, to a fluid of icosahedral clusters. 

At high temperatures, the self-assembling van der Waals 
fluid exhibits liquid-vapour coexistence behaviour exactly 
analogous to the pure van der Waals fluid, and the binodals 
of the self-assembling fluid follow those of the van der Waals 



fluid. This liquid-vapour phase separation is driven by the at- 
tractions between monomers. In this region, the phase bound- 
ary exhibits a positive gradient in the p-T phase diagram 
(Fig. 6). The Clapeyron equation relates this gradient to 



dp A trs H AfrsS 



(19) 



and since the entropy, enthalpy and volume changes are all 
negative for the vapour-liquid transition, the resulting slope is 
naturally positive. 

By contrast, at sufficiently low temperatures, T c \ us (p) lies 
at lower density than the van der Waals binodal. Therefore, 
as the density increases, the vapour forms clusters before it 
reaches the point at which it would otherwise have become 
phase separated. As clusters have no attractions between 
them, the vapour is now stable and does not phase separate. 

As the density of this system approaches that for the max- 
imum density of clusters, packing constraints become impor- 
tant. If the density is to move beyond 12p c [ us max and re- 
main homogeneous, the system must form a monomer-cluster 
mixture with the proportion of monomers increasing with p. 
However, unlike in subsection II C, the system now prefers to 
demix into a vapour of clusters and a monomeric liquid be- 
cause of a lack of attractions between clusters and monomers. 
Consequently, both the solubility of monomers in the cluster 
fluid and of clusters in the monomeric liquid is very low. For 
example, at kT/e = 0.08, p mon cr 3 = 3.44 x 10~ 7 at the lower- 
density binodal (pa 3 = 0.498), and p dus cr 3 = 7.28 x 10~ 120 
at the higher-density binodal (pa 3 = 0.955). This demix- 
ing is a density-driven transition, and in this regime the po- 
sition of the lower-density binodal varies little with temper- 
ature. As the vapour phase is now lower in energy than the 
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FIG. 6. The p-T phase diagram for the self-assembling van der Waals 
fluid at different values of B/b. The critical point is labelled 'CP'. 
The clustering temperature (T c i us ) is also plotted for B/b = 24 (thick 
dottedline) in the monomer-cluster region only, a/ea 3 = l,b/a 3 = 



1, a„ 



0.5. 



liquid due to the strong internal bonding in the clusters, the en- 
thalpy change associated with vapour-liquid transition is now 
positive, while the volume change is of course still negative. 
Hence, the cluster-liquid transition has a negative slope in the 
p-T diagram (Fig. 6). 

The crossover between these two regimes occurs at temper- 
atures near the point at which T c i us (p) crosses the pure van der 
Waals binodal and leads to a rapid increase in the density of 
the lower-density binodal with decreasing temperature until 
12p c ius, max is approached. The width of this crossover regime 
reflects the width of the monomer-cluster equilibrium. The 
change in slope of the binodal occurs slightly above T c i us , as 
the vapour first begins to cluster. Initially, when monomers 
are still in the majority at the binodal, the transition is still 
driven by the attractions between the monomers, and the bin- 
odal occurs when the monomer density reaches the density of 
the binodal for the pure van der Waals system. However, as 
clusters begin to predominate in the coexisting mixture, the 
phase transition is driven more by the lack of attractions be- 
tween monomers and clusters, and occurs when the monomer 
density reaches a critical solubility. 

These changes in the phase diagram are also reflected in the 
p-v and p-p diagrams (Fig. 4). For 0.23 e/k ;$ T < T cr j t , 
we observe van der Waals-like p-v and p-p curves. How- 
ever, below this temperature, even though not yet in the region 
of reentrance, we can observe additional metastable closed 
curves in both p-v and p-p plots. As the temperature de- 
creases, the cluster-rich loop (as seen at kT/e = 0.21) merges 
with the van der Waals curve and thus exhibits two addi- 
tional spinodal points compared to the situation in the pure 
van der Waals case, due to a back-bending of the p-v curve. 
This merging occurs when T c i us (p) meets the spinodal den- 
sity, as shown in Fig. 5. Below this temperature, because the 
cluster predominates in the vapour, the vapour is no longer 
destabilised by the presence of attractions and instead only 
becomes unstable when the density goes beyond 12p c i us max . 
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FIG. 7. Binodal lines for (a) B/b = 12.5 and (b) B/b = 12 are 
shown in violet, compared with the dotted red line representing the 
pure van der Waals binodal. The labels 'c', 'm', T and T refer 
to cluster gas, monomer gas, liquid and fluid states, respectively. 
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To obtain spinodals, we numerically calculated derivatives 
obtained directly from p-v curves. As already mentioned 
in subsection III C, the unusual sections of backbending in 
the p-v isotherms (e.g. curve BC in Fig. 4(c)) are, although 
mechanically stable, compositionally unstable, and therefore 
the points at which dp/dv = oo are not 'true' borders of 
(in)stability. The physically relevant spinodals are those of 
type A and B in the van der Waals-like region of the phase 
diagram, and then D instead of B in lower regions; these are 
shown in Fig. 5. 

For B/b — 24, despite the reentrance, at all temperatures 
below T cr i t there is always a density range for which phase 
separation occurs. We will now consider how the topology of 
the phase diagram can be modified by varying the parameters 
of the model. In particular, we want to investigate whether 
the phase separation at low temperature can be removed, and 
if so, whether the phase diagram might exhibit a closed-loop 
similar to the nicotine-water mixture. 

First, we investigate the effect of decreasing B, as this pa- 
rameter determines p c \ u ^ max and hence the onset of demixing 
in the low-temperature regime. As B/b decreases, the clus- 
ters are able to pack more efficiently and the lower-density 
binodal moves to higher density. Initially, this does not lead 
to any change in the topology of the phase diagram. How- 
ever, the situation becomes more interesting as the binodal 
associated with the cluster fluid meets and crosses that as- 
sociated with the monomelic liquid. Fig. 7(a) illustrates the 
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type of phase diagram that results. We note that this scenario 
does not lead to a lower critical point; rather, at point A in the 
phase diagram, there are two coexisting fluids with the same 
density, but which differ significantly in other characteristics. 
One fluid is a monomeric liquid with high entropy, while the 
other is a dense cluster fluid with low energy. These two flu- 
ids do not evolve continuously into each other as point A is 
approached (as would need to be the case for there to be a 
critical point), but instead retain their separate identities. 

As AV = at point A in Fig. 7, dp/ 'dT — !• oo and there 
is a change in sign of the slope of the phase boundary in the 
p-T phase diagram. Initially, beyond this point, the coexisting 
cluster fluid is more dense, and so dp /dT is again positive. 
However, for B/b > 12, the clusters still pack less efficiently, 
and so at low temperature and densities beyond 12p c i USjmax , 
there must again be a region of demixing because of the lack 
of attractions between cluster and monomers. The slope of 
the p-T phase boundary corresponding to this demixing must 
again be negative and so there must be a second point (B) 
where the slope of the phase-boundary changes and where the 
coexisting fluids have the same density. Curiously, for this 
phase diagram topology, there is a temperature range (0.19 
kT/e < 0.21 for B/b = 12.5) where there are three fluicU 
fluid phase transitions as the system is pressurised. 

Finally, at B/b = 12, there must be a further change in 
the phase diagram topology, because at this point monomers 
no longer pack more efficiently than clusters, and at low tem- 
perature the cluster fluid will be more stable for any feasible 
density. As B/b approaches 12 from above, point B moves 
to higher and higher temperatures, and the density range for 
the cluster-fluid/monomeric liquid demixing becomes smaller 
and smaller until it disappears at B/b = 12. The resulting 
phase diagram is illustrated in Fig. 7(b). 

The second effect we wish to consider is changing the po- 
sition of the liquid-vapour critical point with respect to the 
clustering transition. For a given b, T cl -i t is determined by the 
value of a, and T c i us is determined by cr pw . In contrast to our 
model, for the simulated patchy particle systems, the two pa- 
rameters are not independent, but instead the position of the 
critical point varies with the patch width (and the number of 
patches). 2S,:,,S For example, it can be seen from Fig. 1 that the 
temperature range associated with liquid-vapour coexistence 
decreases as cr pw decreases. Here, we therefore choose to vary 
one of the parameters (a) whilst keeping the other one (cr pw ) 
fixed. 

The results for a variety of values of a are shown in 
Fig. 8. When the van der Waals critical temperature is above 
the clustering temperature (for example, at a = 1 ea 3 and 
a = 0.8 £<t 3 , which correspond to critical temperatures of 
kTcnt/e = 0.296 and kT aa /e — 0.237, respectively), van der 
Waals behaviour is observed at high temperatures. When the 
critical temperature descends below the clustering tempera- 
ture, we nevertheless continue to observe a bending of the bin- 
odal curve to lower densities in the region where liquid-vapour 
coexistence previously took place. However, this feature is not 
some remnant of the liquid-vapour transition, but is instead 
associated with demixing in the vicinity of T c \ us (p). At T c i us , 
there is an equal probability of particles being monomeric or 
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FIG. 8. A T-p phase diagram for the self-assembling van der Waals 
fluid with a selection of values for the van der Waals attractive pa- 
rameter a. b/a 3 = 1, b/B = 24, a pw = 0.5. 



as part of clusters; however, if p and a are sufficiently large, 
the system will not form a homogeneous mixture, but will 
instead demix. The narrowness in temperature of this fea- 
ture is associated with the sharpness with which the chemi- 
cal equilibrium shifts between being monomer dominated to 
being cluster dominated as the temperature decreases. As a 
decreases, there is a progressively smaller demixing driving 
force and so the size of this feature decreases. Note that these 
changes lead to a non-monotonic dependence of T cl i t on a, 
because T a -[ t initially increases as p a \ t increases, mirroring the 
dependence of T c i us on p (Fig. 2). As a becomes smaller still, 
the clusters and monomers begin to be able to mix at high 
temperatures because less entropy of mixing is required to 
overcome the enthalpic mixing incompatibility. Hence, T cr i t 
again begins to decrease. However, even at a small value of a, 
e.g. a = 0.05 ecr 3 in Fig. 8, there is still a substantial region 
of demixing. As a tends to zero, the demixing regime disap- 
pears, recovering the monomer-cluster equilibrium studied in 
section II. 

So far, we have exclusively considered particles that can 
self-assemble into icosahedral clusters, but our approach can 
easily be generalised to the formation of other types of clus- 
ter. Here, we illustrate this by calculating the phase diagram 
for particles able to form tetrahedra. The requisite changes 
in the model from the icosahedral system discussed above are 
as follows. The moment of inertia of a tetrahedron about the 
principal axes is / = 4m<r 2 , and the mass of the cluster is 
Am. The symmetry number of a tetrahedron is 12, and that 
of an individual monomer is 3. The ground state energy of 
the tetrahedron is — 6e. Finally, the circumradius of a tetrahe- 
dron with edge length a is r = ^/3/8 a, and so the phys- 
ical value of B/b which we have used is B/b = 11. This 
gives a maximum particle density associated with clusters of 
4p c ius, maxC 3 = 4/11 s» 0.36. We have assumed that tetrahe- 
dra will pack approximately like spheres; this approximation 
is not so far-fetched as it might seem at first glance, and recent 
research confirms that the packing fractions are very similar. 
We also scale the value of the parameter a by 3/5 compared 
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FIG. 9. A T-p phase diagram for the self-assembling van der Waals 
fluid forming tetrahedral clusters. The solid line is the binodal curve, 
whilst the thin dotted line represents the pure van der Waals system 
binodal. The thick dashed line is the clustering temperature, T c i m (p), 
for this system. The labels 'c', 'm', '1' and 'f refer to cluster gas, 
monomer gas, liquid and fluid states, respectively, a/ea 3 = 0.6, 
B/b = 11, b/a 3 = 1, cr pw = 0.5. 

to that for the icosahedron in order to take into account that 
there are only three rather than five patches that give rise to 
the interactions that drive the liquid- vapour phase separation. 
The phase diagram of this system is depicted in Fig. 9. The 
behaviour of the system is similar to the icosahedral system, 
with the main differences being easy to rationalise. For exam- 
ple, the clustering temperature is significantly lower due to the 
lower energetic driving force to form clusters. Similarly, the 
crossover between the liquid-vapour and demixing regimes 
occurs over a wider temperature range due to the greater width 
of the monomer-cluster equilibrium due to the smaller size of 
the clusters. 



IV. DISCUSSION AND CONCLUSIONS 

We have calculated the complete phase diagram for a sys- 
tem exhibiting competition between phase separation and self- 
assembly into monodisperse clusters. For nearly all of the 
parameterisations of the model that we have considered, de- 
creasing the temperature at constant density results, for a lim- 
ited range of density, in the transition from a gas of monomers, 
to a liquid-gas phase-separated mixture, and then to a gas 
of clusters. Therefore, the existence of reentrant phase be- 
haviour, driven by the lower energy of clusters, is clear within 
this model and is a robust feature of such systems. In this as- 
pect, the model agrees with our group's simulation results for 
patchy colloidal particles, • and the Janus particle simula- 
tions of Sciortino and co-workers, 33 ' 37 both of which showed 
evidence of reentrant phase behaviour. 

In particular, a comparison of the final phase diagram with 
the yield plot shown in Fig. 1 shows that decreasing the tem- 
perature at a constant density of per 3 = 0.15 and a patch 
width of <7 pw = 0.5 involves a reentrant phase transition in 
both cases. Even though we have not chosen to match the pa- 
rameter a precisely to the icosahedral system studied in sim- 
ulations, there is surprisingly good agreement not only in the 
underlying physical behaviour, but also in the numerical val- 



ues. 

While the phase behaviour of the fluid at low temperature 
is not yet fully clear in the simulated systems, our model al- 
lows the complete phase diagram to be calculated. We have 
established that phase separation persists at low temperatures 
within the framework of our model. The driving force for it, 
however, is not the attractions between monomers as in the 
pure van der Waals fluid, but rather cluster-liquid immiscibil- 
ity. For this reason, reentrance is observed over only a limited 
range of density in all physically relevant phase diagrams; at 
higher densities, demixing between the cluster and the liquid 
states ensures continued phase separation. 

In our model, the change in topology of the phase diagram 
as the model parameters are varied is easy to calculate. The 
free volume parameter ratio B/b changes the maximum den- 
sity of the cluster, and therefore affects the position of the 
lower-density binodal of the T-p phase diagram in the demix- 
ing regime. The change in the van der Waals parameter a, on 
the other hand, changes the van der Waals critical temperature, 
and, provided that it is greater than the clustering temperature, 
affects the temperature range of the phase-separated region. 

Of course, in our theoretical model, there is a significant 
number of approximations, the effect of which we should 
try to understand, in particular in relation to the compari- 
son of our results with phase diagrams obtained from simu- 
lations. Firstly, our focus has just been on the fluid behaviour 
of these self-assembling systems, and, unlike in simulations, 
we can explore the low temperature form of the fluid phase 
behaviour without having to worry about other phases obscur- 
ing the fundamental behaviour. However, crystallisation (of 
both monomers and clusters) is likely to have a major effect 
on the overall form of the phase diagram. There are prece- 
dents for such cluster crystals; for example, icosahedral virus 
particles (monodisperse aggregates of virus capsid proteins, 
perhaps with a nucleic acid genome inside) are often observed 
to form crystals when they are densely packed in cells. 

We can construct a schematic phase diagram illustrating 
the potential effects of crystallisation on it by assuming hard- 
sphere crystallisation of both components. This approxima- 
tion is probably quite reasonable for the clusters, as all the 
attractive patches are involved in the internal bonding of the 
cluster, but probably less so for the monomers, where the ef- 
fect of the patches on the crystal stability and form is less 
clear. For hard spheres, the fluid-crystal coexistence limits 
occur at packing fractions <f> = 0.494 and </> = 0.545, 61 and 
by associating the maximum densities of the fluids with the 
packing fraction in random close packing, </> = 0.64, 62 we 
can obtain estimates for the densities of the fluid-crystal bin- 
odals. Further, the maximum hard-sphere crystalline packing 
fraction, cj) = 0.74, provides a means for obtaining an upper 
density limit for the stability of the cluster crystal. Finally, the 
clustering temperature provides an upper temperature stabil- 
ity limit for the cluster crystal. The resulting schematic phase 
diagram is shown in Fig. 10. Note that this phase diagram 
has a second triple point associated with the disappearance of 
the cluster crystal phase, as well as the usual one associated 
with the disappearance of the liquid phase. We have not at- 
tempted to consider the effect of other possible phases on the 
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FIG. 10. A schematic T-p phase diagram to illustrate the potential 
effects of crystallisation on the phase diagram in Fig. 5. 



phase diagram, such as the lamellar phase observed for Janus 
particles 33 ' 37 and the two-dimensional crystalline lamellae 
seen in simulations of octahedron-forming patchy particles. 2 

Secondly, there are also a series of approximations associ- 
ated with our description of the fluid. For example, the free 
volume approximation is not flexible enough to fully cap- 
ture both the high- and the low-density behaviour (we choose 
B and b values appropriate to the high-density limit) and the 
form of the available volume term differs from that for the 
exact hard-sphere second virial coefficient for interactions be- 
tween monomers and clusters. There is also an obvious prob- 
lem when the density of clusters approaches the maximum 
allowed, as the free volume approximation predicts that the 
translational partition function of both clusters and monomers 
should tend to zero. However, because of the large differ- 
ences in size between monomers and icosahedral clusters, 
even when the clusters are close-packed, there is still space 
available for the monomers in the gaps between the clus- 
ters. In other words, the translational partition function of the 
monomers should not tend to zero as the translational partition 
function of the clusters does so. However, such correlation ef- 
fects associated with the positions of the particles are not so 
straightforward to take into account. 

In our model, we also effectively assume that the clus- 
ters pack, in essence, like hard spheres, without any change 
of shape or compression, regardless of the overall density 
of the system. This assumption, which is responsible for 
the near vertical slope of the lower-density binodal in the 
low-temperature demixing regime, may be reasonable for 
monodisperse clusters that have a well-defined shape and are 
held together by specific bonds (although for our patchy parti- 
cle potential, the clusters will be slightly compressible). How- 
ever, this assumed cluster incompressibility may be a factor 
in the mismatch in the detailed shape of our phase diagrams 
compared to those obtained in the simulations by Sciortino 
and co-workers, • where the shapes of clusters are much 
more deformable, as well as not being monodisperse. In their 
computed phase diagrams, they observe the binodals moving 
to higher densities with decreasing temperature considerably 
more gradually than we do. 

Another approximation we have made is that we can model 
the monomeric liquid and vapour phases as a van der Waals 



fluid. In this model, the critical density has a fixed value (de- 
pending only on the free volume parameter b), whereas we 
know that in patchy particle colloids, the critical point varies 
with both the patch width and the patch number. ■ Never- 
theless, this deficiency is probably not so important for the 
qualitative shape of the phase diagram, especially as we are 
most concerned in this paper with the behaviour at lower tem- 
perature associated with the onset of clustering. 

Finally, we have assumed that clusters form either com- 
pletely or not at all, i.e. only a single size of clusters can be 
formed. Although this is an approximation, the probability of 
observing partly-formed clusters at equilibrium is very small, 
and so the effect on the phase diagrams will be negligible. 

Although these approximations will affect the quantitative 
comparison between the phase diagrams for the current model 
and those computed in patchy particle simulations, 29 ' 32 and 
hopefully at some future point those found experimentally for 
self-assembled patchy colloids, we believe that the fundamen- 
tal insights the current model gives will help to provide a theo- 
retical framework for understanding the fluid phase behaviour 
of these kinds of system. 
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